rm(list=ls())
devtools::install_github('xuyiqing/fastplm')
devtools::install_github('xuyiqing/fect')
devtools::install_github('xuyiqing/panelview')


################################################################################
#                                    Rondonia ATU
#                                        GFC
#                                     BZ 20 km
################################################################################

### prepare data

panel_GFC <- read_dta('.../DIDmatching_BZ20km_GFC.dta')

# subset of data for rondonia
subset <- panel_GFC[,c("ID1", "year","enacR_10","agrifires" ,"state_", "tavg", "light", "prec","treecoverloss", "ID_state", "rdo10PSM_ron_weight", "rdo10PSM_rons_weight", "rondonia_weight", "PA0108")]

# subset of untreated group and matched treated
rons_ATU <- subset[is.na(subset$PA0108) & subset$ID_state==22 & (subset$enacR_10==0 | subset$enacR_10==1 & !is.na(subset$rondonia_weight)),]
rons_ATU<- rons_ATU[!is.na(rons_ATU$ID1),]
rons_ATU <- rons_ATU[,c("ID1", "year","enacR_10","agrifires" ,"state_", "tavg", "light", "prec","treecoverloss", "ID_state", "rondonia_weight" )]

# subset for state size reduction
rons_ATU<-rons_ATU %>% group_by(ID1) %>% mutate(enacR_10_ = max(enacR_10))
rons_ATU<-rons_ATU[(rons_ATU$enacR_10_==1 & rons_ATU$state_==1) | (rons_ATU$enacR_10==0 & is.na(rons_ATU$state_)) ,]
rons_ATU<- rons_ATU[!is.na(rons_ATU$ID1),]

### DID estimation and plot
fect_rons <- fect(treecoverloss ~ enacR_10 + agrifires + tavg+ light + prec, data =rons_ATU, 
                   index = c("ID1","year"), 
                   method = "fe", force = "two-way", se = TRUE)

plot(fect_rons, main = "Estimated ATU", ylab = "Effect of PA Size Reduction on Tree Cover Losses", 
     cex.main = 0.8, cex.lab = 0.8, cex.axis = 0.8)

### test of parallel trend 
plot(fect_rons, type = "equiv", ylim = c(-0.05,0.05), 
     cex.legend = 0.6, main = "Testing Pre-Trend (FEct)", cex.text = 0.8, pre.periods = c(-4,0))


################################################################################
#                                    Rondonia ATU
#                                        GFC
#                                     BZ 10 km
################################################################################

panel_GFC_bz10 <- read_dta('.../DIDmatching_BZ10km_GFC.dta')

### prepare data

# subset of data for rondonia
subset <- panel_GFC_bz10[,c("ID1", "year","enacR_10_B3","agrifires" ,"state_", "tavg", "light", "prec","treecoverloss", "ID_state", "rondonia_weight", "PA0108")]

# subset of untreated group and matched treated
rons_ATU_bz10 <- subset[is.na(subset$PA0108) & subset$ID_state==22 & (subset$enacR_10_B3==0 | subset$enacR_10_B3==1 & !is.na(subset$rondonia_weight)),]
rons_ATU_bz10<- rons_ATU_bz10[!is.na(rons_ATU_bz10$ID1),]
rons_ATU_bz10 <- rons_ATU_bz10[,c("ID1", "year","enacR_10_B3","agrifires" ,"state_", "tavg", "light", "prec","treecoverloss", "ID_state", "rondonia_weight" )]

# subset for state size reduction
rons_ATU_bz10<-rons_ATU_bz10 %>% group_by(ID1) %>% mutate(enacR_10_B3_ = max(enacR_10_B3))
rons_ATU_bz10<-rons_ATU_bz10[(rons_ATU_bz10$enacR_10_B3_==1 & rons_ATU_bz10$state_==1) | (rons_ATU_bz10$enacR_10_B3==0 & is.na(rons_ATU_bz10$state_)) ,]
rons_ATU_bz10<- rons_ATU_bz10[!is.na(rons_ATU_bz10$ID1),]

### DID estimation 
fect_rons_bz10 <- fect(treecoverloss ~ enacR_10_B3 + agrifires + tavg+ light + prec, data =rons_ATU_bz10, 
                   index = c("ID1","year"), 
                   method = "fe", force = "two-way", se = TRUE)

### test of parallel trend 
plot(fect_rons_bz10, type = "equiv", ylim = c(-0.05,0.05), 
     cex.legend = 0.6, main = "Testing Pre-Trend (FEct)", cex.text = 0.8, pre.periods = c(-4,0))


################################################################################
#                                    Rondonia ATU
#                                        TMF
#                                     BZ 20 km
################################################################################


panel_TMF <- read_dta('.../DIDmatching_BZ20km_TMF.dta')

### prepare data

# subset of data for rondonia

subset <- panel_TMF[,c("ID1", "year","enacR_10","agrifires" ,"state_", "tavg", "light", "prec","TMF_", "ID_state", "rondonia_weight", "PA0108")]

# subset of untreated group and matched treated

rons_ATU_tmf <- subset[is.na(subset$PA0108) & subset$ID_state==22 & (subset$enacR_10==0 | subset$enacR_10==1 & !is.na(subset$rondonia_weight)),]
rons_ATU_tmf<- rons_ATU_tmf[!is.na(rons_ATU_tmf$ID1),]
rons_ATU_tmf <- rons_ATU_tmf[,c("ID1", "year","enacR_10","agrifires" ,"state_", "tavg", "light", "prec","TMF_", "ID_state", "rondonia_weight" )]

# subset for state size reduction

rons_ATU_tmf<-rons_ATU_tmf %>% group_by(ID1) %>% mutate(enacR_10_ = max(enacR_10))
rons_ATU_tmf<-rons_ATU_tmf[(rons_ATU_tmf$enacR_10_==1 & rons_ATU_tmf$state_==1) | (rons_ATU_tmf$enacR_10==0 & is.na(rons_ATU_tmf$state_)) ,]
rons_ATU_tmf<- rons_ATU_tmf[!is.na(rons_ATU_tmf$ID1),]

### DID estimation 

fect_rons_tmf <- fect(TMF_ ~ enacR_10 + agrifires + tavg+ light + prec, data =rons_ATU_tmf, 
                       index = c("ID1","year"), 
                       method = "fe", force = "two-way", se = TRUE)


### test of parallel trend 

plot(fect_rons_tmf, main = "Estimated ATU", ylab = "Effect of PA Size Reduction on Deforestation and Forest Degradation", 
     cex.main = 0.8, cex.lab = 0.8, cex.axis = 0.8)
print(fect_rons_tmf)



################################################################################
#                                    Para ATT
#                                        GFC
#                                     BZ 20 km
################################################################################

### prepare data

# subset of data for Para

subsetpara <- panel_GFC[,c("ID1", "year","enacR_12","agrifires" , "tavg", "light", "prec","treecoverloss", "ID_state", "para_weight", "PA0108")]

# subset of treated group and matched untreated

para_ATT <- subsetpara[is.na(subsetpara$PA0108) & subsetpara$ID_state==14 & (subsetpara$enacR_12==1 | subsetpara$enacR_12==0 & !is.na(subsetpara$para_weight)),]
para_ATT<- para_ATT[!is.na(para_ATT$ID1),]
para_ATT <- para_ATT[,c("ID1", "year","enacR_12","agrifires" , "tavg", "light", "prec","treecoverloss",  "para_weight" )]

### DID estimation 

fect_para <- fect(treecoverloss ~ enacR_12 + agrifires + tavg+ light + prec, data =para_ATT, 
                   index = c("ID1","year"), 
                   method = "fe", force = "two-way", se = TRUE)

### test of parallel trend 

plot(fect_para, type = "equiv", ylim = c(-0.05,0.05), 
     cex.legend = 0.6, main = "Testing Pre-Trend (FEct)", cex.text = 0.8)


################################################################################
#                                    Para ATT
#                                        GFC
#                                     BZ 10 km
################################################################################

### prepare data

# subset of data for Para

subsetpara <- panel_GFC_bz10[,c("ID1", "year","enacR_12_B3","agrifires" , "tavg", "light", "prec","treecoverloss", "ID_state", "para_weight", "PA0108")]

# subset of treated group and matched untreated

para_att_bz10 <- subsetpara[is.na(subsetpara$PA0108) & subsetpara$ID_state==14 & (subsetpara$enacR_12_B3==1 | subsetpara$enacR_12_B3==0 & !is.na(subsetpara$para_weight)),]
para_att_bz10<- para_att_bz10[!is.na(para_att_bz10$ID1),]
para_att_bz10 <- para_att_bz10[,c("ID1", "year","enacR_12_B3","agrifires" , "tavg", "light", "prec","treecoverloss",  "para_weight" )]

### DID estimation 

fect_para_bz10 <- fect(treecoverloss ~ enacR_12_B3 + agrifires + tavg+ light + prec, data =para_att_bz10, 
                   index = c("ID1","year"), 
                   method = "fe", force = "two-way", se = TRUE)

### test of parallel trend 

plot(fect_para_bz10, type = "equiv", ylim = c(-0.05,0.05), 
     cex.legend = 0.6, main = "Testing Pre-Trend (FEct)", cex.text = 0.8)


################################################################################
#                                    Para ATT
#                                        TMF
#                                     BZ 20 km
################################################################################

### prepare data

# subset of data for Para

subsetpara <- panel_TMF[,c("ID1", "year","enacR_12","agrifires" , "tavg", "light", "prec","TMF_", "ID_state", "para_weight", "PA0108")]

# subset of treated group and matched untreated

para_att_tmf <- subsetpara[is.na(subsetpara$PA0108) & subsetpara$ID_state==14 & (subsetpara$enacR_12==1 | subsetpara$enacR_12==0 & !is.na(subsetpara$para_weight)),]
para_att_tmf<- para_att_tmf[!is.na(para_att_tmf$ID1),]
para_att_tmf <- para_att_tmf[,c("ID1", "year","enacR_12","agrifires" , "tavg", "light", "prec","TMF_", "ID_state", "para_weight" )]

### DID estimation

fect_para_tmf <- fect(TMF_ ~ enacR_12 + agrifires + tavg+ light + prec, data =para_att_tmf, 
                       index = c("ID1","year"), 
                       method = "fe", force = "two-way", se = TRUE)

### test of parallel trend 

plot(fect_para_tmf, type = "equiv", ylim = c(-0.05,0.05), 
     cex.legend = 0.6, main = "Testing Pre-Trend (FEct)", cex.text = 0.8, pre.periods = c(-4,0))


################################################################################
#                                   Roraima ATU
#                                        GFC
#                                     BZ 20 km
################################################################################

### prepare data

# subset of data for Roraima

subsetror <- panel_GFC[,c("ID1", "year","enacR_09","agrifires" , "tavg", "light", "prec","treecoverloss", "ID_state", "roraima_weight", "PA0108")]

# subset of untreated group and matched treated

ror_atu_gfc <- subsetror[is.na(subsetror$PA0108) & subsetror$ID_state==23 & (subsetror$enacR_09==0 | subsetror$enacR_09==1 & !is.na(subsetror$roraima_weight)),]
ror_atu_gfc<- ror_atu_gfc[!is.na(ror_atu_gfc$ID1),]
ror_atu_gfc <- ror_atu_gfc[,c("ID1", "year","enacR_09","agrifires" , "tavg", "light", "prec","treecoverloss",  "roraima_weight" )]

### DID estimation 

fect_ror <- fect_ror(treecoverloss ~ enacR_09 + agrifires + tavg+ prec, data =ror_atu_gfc, 
                     index = c("ID1","year"), 
                     method = "fe", force = "two-way", se = TRUE)

### test of parallel trend 

plot(fect_ror, type = "equiv", ylim = c(-0.05,0.05), 
     cex.legend = 0.6, main = "Testing Pre-Trend (fect_ror_gfc)", cex.text = 0.8)


################################################################################
#                                   Roraima ATU
#                                        GFC
#                                     BZ 10 km
################################################################################

### prepare data

# subset of data for Roraima

subsetror <- panel_GFC_bz10[,c("ID1", "year","enacR_09_B3","agrifires" , "tavg", "light", "prec","treecoverloss", "ID_state", "roraima_weight", "PA0108")]

# subset of untreated group and matched treated

ror_atu_bz10 <- subsetror[is.na(subsetror$PA0108) & subsetror$ID_state==23 & (subsetror$enacR_09_B3==0 | subsetror$enacR_09_B3==1 & !is.na(subsetror$roraima_weight)),]
ror_atu_bz10<- ror_atu_bz10[!is.na(ror_atu_bz10$ID1),]
ror_atu_bz10 <- ror_atu_bz10[,c("ID1", "year","enacR_09_B3","agrifires" , "tavg", "light", "prec","treecoverloss",  "roraima_weight" )]

### DID estimation 

fect_ror_bz10 <- fect(treecoverloss ~ enacR_09_B3 + agrifires + tavg+ prec, data =ror_atu_bz10, 
                     index = c("ID1","year"), 
                     method = "fe", force = "two-way", se = TRUE)


### test of parallel trend 

plot(fect_ror_bz10, type = "equiv", ylim = c(-0.05,0.05), 
     cex.legend = 0.6, main = "Testing Pre-Trend (FEct)", cex.text = 0.8)

################################################################################
#                                   Roraima ATU
#                                        TMF
#                                     BZ 20 km
################################################################################

### prepare data

# subset of data for Roraima

subsetror <- panel_tmf[,c("ID1", "year","enacR_09","agrifires" , "tavg", "light", "prec","TMF_", "ID_state", "roraima_weight", "PA0108")]

# subset of untreated group and matched treated

ror_atu_tmf <- subsetror[is.na(subsetror$PA0108) & subsetror$ID_state==23 & (subsetror$enacR_09==0 | subsetror$enacR_09==1 & !is.na(subsetror$roraima_weight)),]
ror_atu_tmf<- ror_atu_tmf[!is.na(ror_atu_tmf$ID1),]
ror_atu_tmf <- ror_atu_tmf[,c("ID1", "year","enacR_09","agrifires" , "tavg", "light", "prec","TMF_",  "roraima_weight" )]

### DID estimation and plot

fect_ror_tmf <- fect(TMF_ ~ enacR_09 + tavg+ prec, data =ror_atu_tmf, 
                         index = c("ID1","year"), 
                         method = "fe", force = "two-way", se = TRUE)


### test of parallel trend 

plot(fect_ror_tmf, type = "equiv", ylim = c(-0.05,0.05), 
     cex.legend = 0.6, main = "Testing Pre-Trend (FEct)", cex.text = 0.8)


